Three-dimensional nonlinear model of rock creep under freeze–thaw cycles

In areas with large differences between day and night temperature, the freeze–thaw cycle and frost heaving force in rock mass generate cracks within the rock, which seriously threatens the stability and safety of geotechnical engineering structures and surrounding buildings. This problem can be solved by developing a reasonable model that accurately represents the rock creep behavior. In this study, we developed a nonlinear viscoelastic–plastic creep damage model by introducing material parameters and a damage factor while connecting an elastomer, a viscosity elastomer, a Kelvin element, and a viscoelastic–plastic element in series. One- and three-dimensional creep equations were derived, and triaxial creep data were used to determine the model parameters and to validate the model. The results showed that the nonlinear viscoelastic–plastic creep damage model can accurately describe rock deformation in three creep stages under freeze–thaw cycles. In addition, the model can describe the time-dependent strain in the third stage. Parameters G1, G2, and η20’ decrease exponentially with the increase in the number of freeze–thaw cycles while parameter λ increases exponentially. These results provide a theoretical basis for studying the deformation behavior and long-term stability of geotechnical engineering structures in areas with large diurnal temperature differences.


Introduction
In the past few years, the number and scale of open-pit mines in China have increased rapidly. The number of open-pit mines in Inner Mongolia, Xinjiang, Gansu, and other regions in Northwest China was 83.39% of the national count in 2021. Rock mass engineering in these areas is faced with complex geographical and climatic conditions such as seasonally frozen soil and low temperatures [1]. The large variations in temperature due to seasonal alternations and the day-night cycle result in varying degrees of freeze-thaw disasters such as frost heave cracking, freeze-thaw collapse, and slope instability [2,3]. These disasters directly affect rock engineering stability and safety and pose serious challenges to engineering construction and subsequent maintenance [4][5][6]. Therefore, in-depth studies on the aging mechanical behavior and creep law of rock specimens under freeze-thaw environments are of vital theoretical significance for enriching the creep mechanics theory and analyzing the time-dependent strain of rocks. Previous research has mostly focused on the physical and mechanical behavior and creep deformation law of rocks under freeze-thaw cycles [7][8][9], and important results have been obtained through experimental analyses, quantitative evaluations, and engineering practice. According to the Lemaitre strain equivalence hypothesis, Yuan et al. [10] set the initial damage condition of intact rock as the reference damage condition to construct a damage constitutive model. Yan et al. [11] separated the strain of rock under the action of freeze-thaw cycles into initial damage strain, additional damage strain, and plastic strain, and the authors established a damage constitutive model of the elastic-plastic freeze-thaw cycle. Considering the theory of continuous damage mechanics, Zhang et al. [12] constructed a damage model of freeze-thaw loading, factoring in changes in confining pressure. In terms of the Lemaitre strain equivalence hypothesis, Chen et al. [13] established a jointed rock mass composite damage model, considering macroscopic and microscopic defects. Liu et al. [14] used the definition of damage and the Lemaitre strain equivalent equation to establish a freeze-thaw damage evolution equation in the constant creep stage. According to the von Mises yield criterion, a damage formula for freeze-thaw and creep changes throughout the nonlinear creep stage was constructed. The Burgers model parameters were corrected by introducing damage variables to construct the creep damage model under the action of freeze-thaw cycles. Li [15] constructed a rock freezethaw damage element based on nuclear magnetic resonance porosity and a creep damage element based on volumetric strain; both elements were introduced into classical rheological elements to construct a nonlinear creep damage model for describing the freeze-thaw behavior of rocks. Zhang et al. [16] conducted shear creep experiments and mesoscopic characteristics analyses, proposed rock freeze-thaw damage viscous components, and established a granite freeze-thaw shear creep constitutive model for characterizing the influence of freeze-thaw cycles on shear creep behavior. Through a triaxial creep test, Wan et al. [17] established a Nishihara model based on damage mechanics principles and extended the model to a threedimensional form according to the elastic-plastic mechanics theory. Song et al. [18] performed triaxial creep experiments on saturated red sandstone under freeze-thaw cycles and established a freeze-thaw damage creep model considering both freeze-thaw cycles and damage effects. Hou et al. [19] combined statistical damage mechanics theory and the maximum tensile strain failure criterion to establish a freeze-thaw rock damage model. Zhang et al. [20] performed mesoscopic characteristics analyses and shear creep tests on granite specimens undergoing different freeze-thaw cycles. They found that the increase in the number of freeze-thaw cycles increased the damage degree of the rock surface, and a creep model of freeze-thaw coupling damage and stress damage was constructed to describe the rock creep characteristics.
Evidently, key results have been obtained in the study of the mechanical characteristics and creep deformation behavior of rocks under freeze-thaw cycles. In this study, we constructed a nonlinear viscoelastic-plastic creep damage model. Material parameters and the damage factor were introduced into the model, as well as an elastomer, a viscosity elastomer, a Kelvin element, and a viscoelastic-plastic element connected in series. The one-and three-dimensional creep equations were derived, and triaxial creep data were used to obtain the model parameters and validate the model. The findings of this study can provide a scientific basis and technical reference for the support maintenance and long-term stability analysis of geotechnical engineering structures in Northwest China, which has seasonally frozen soil and an extremely lowtemperature complex climate environment.

Analysis of triaxial creep test of mudstone under freeze-thaw cycles
Li [21] conducted triaxial creep experiments of mudstone specimens under 0, 10, 20, and 30 freeze-thaw cycles. In this study, we put all the mudstone specimens into the vacuum saturation equipment for 48 h under negative pressure saturation before conducting the freeze-thaw test. The maximum and minimum temperature limits were 20˚C and -20˚C, respectively. Chen's loading method was used for the triaxial rock creep test, and an LVDT displacement meter was used for displacement monitoring. The confining pressure was maintained at 5 MPa, and the axial loading was divided into five levels: 10 MPa, 20 MPa, 30 MPa, 40 MPa, and 50 MPa. The loading rate was controlled by the force control mode to 30 kN/min approximately 0.25 MPa/s), and the load was applied to each creep stability stage for 18 h.   Under each loading stage, the rock specimens underwent elastic deformation or viscoelastic deformation due to the sudden increase in axial load, which is approximately expressed as a vertical or high slope line segment on the creep curve. The first deformation stage is characterized by a slow increase in mudstone deformation and a decrease in creep deformation. Thus prior to the second stage, the mudstone deformation increases gradually while the creep rate tends to be stable. The first and second stages are characterized by a few internal mudstone cracks and the growth and germination of micro cracks. In these stages, the mudstone damage is minimal, and no obvious damage crack appears on the mudstone surface. When the stress on the rock exceeds the yield strength of the rock, the rock specimen enters the third deformation stage. In this stage, the axial deformation and deformation rate of the mudstone increase rapidly, and macroscopic cracks appear until failure.

Establishment of nonlinear viscoelastic-plastic creep damage model and derivation of creep equation
In the first and second stages, zero or negligible damage occurs in the rock specimen; in the third stage, damage gradually accumulates in the rock and breaks the rock, showing obvious nonlinear characteristics [23,24]. The classical Burgers model and the Nishihara model, which are composed of ideal elements such as an elastomer, a Kelvin element, and a viscosity element, can only model the rock creep behavior rather than the entire creep process. Thus, describing the creep properties of the third stage is a key scientific problem [25]. Therefore, according to existing results of the Burgers model, we modeled the creep behavior in the third stage by considering the time correlation of the viscosity coefficient of the Kelvin element; we also used the viscoelastic-plastic element considering damage (Fig 2).
By summing the strains of each element, the Burgers model [26] creep equation can be expressed as: where σ is the stress, MPa; k 1 is the elastic modulus of the elastomer, MPa; η 1 is the viscosity coefficient of the viscosity element, MPa�h; k 2 is the elastic modulus of the Kelvin element, MPa; and η 2 is the viscosity coefficient of the Kelvin element, MPa�h.
In the first stage, with the increase in time, the unit creep of the rock gradually decreases while the creep curve slope gradually decreases. The entire stage exhibits nonlinear characteristics. In the second stage, the attenuation, characteristics, and time are closely related to the rock stress. During the creep test, the bonding structures in the rock are gradually damaged under constant stress; macroscopically, the creep strain increases gradually. The viscosity parameters increase with time in the first stage but gradually tend to be stable in the second stage, indicating that the rock creep reaches a stable state in the second stage [27,28]. Although the traditional Kelvin element can better depict the creep attenuation law under a given stress, it fails to accurately reflect the nonlinear creep behavior under different stresses because the viscosity parameter remains constant with time. The correlation between the viscosity coefficient and time in the Kelvin element were considered in this study. It is assumed that the coefficient and time in the creep process are power functions [29]. The differential constitutive equation of the nonlinear Kelvin element is ( where η 20 is the initial viscosity coefficient of the nonlinear Kelvin element, MPa�h; η 2 (t) is the viscosity coefficient of the nonlinear Kelvin element, MPa�h 2 ; and λ is a constant. Differentiating Eq (3) yields the creep equation (Eq (4)) of the nonlinear Kelvin element: 3.1.2 Viscoelastic-plastic element considering damage. When the loading level exceeds the long-term length, damage occurs in the rock and gradually develops over time [30]. Compared with the initial and constant creep stages, the accelerated creep stage accounts for a small part of the creep test process. The accelerated creep stage has a short duration but produces significant nonlinear deformation [31]. Therefore, to accurately describe the creep behavior of the third creep stage, we introduce the damage variable to the viscosity coefficient and the elastic modulus to the viscoelastic element to model the rock damage evolution mechanism [32]: Then, ( where D is the damage factor; η 30 is the initial viscosity coefficient of the viscoelastic-plastic element considering damage, MPa�h; η 30 (D,t) is the viscosity coefficient with damage variable, MPa�h; k 30 is the initial elastic modulus of the viscoelastic-plastic element considering damage, MPa; and k 30 (D,t) is the elastic modulus with damage variable, MPa. For simplicity, the classical creep damage function [33] is improved using the simplified damage expression commonly used internationally [34][35][36], and the damage factor D is given as Eq (7) [37,38]: where e is the natural constant, and α is the material constant.
We can obtain the constitutive equation of this element considering damage as Eq (8) [39]: From the integral solution of Eq (8), we can obtain the one-dimensional creep equation of the viscoelastic-plastic element considering damage: Therefore, the one-dimensional creep equation of the nonlinear viscoelastic-plastic creep damage model is: 8 > > > < > > > :

Three-dimensional creep equation
In engineering practice, rock mass is typically in three dimensions with a complex stress state. Hence, the one-dimensional creep equation should be extended to the three-dimensional space to establish a three-dimensional creep equation, providing scientific support for accurately describing rock creep behavior [37]. Geotechnical engineering research has produced several key findings through element models. However, certain problems exist in constructing the three-dimensional rock creep equation and determining parameters of the test curve. These problems are as follows: (1) when the test curve of creep deformation under threedimensional stress conditions is fitted with a one-dimensional constitutive equation, the model parameters have different physical meanings in one-dimensional and three-dimensional cases, and the parameters obtained are inaccurate because a one-dimensional constitutive equation does not consider the effect of confining pressure on rock creep characteristics.
(2) The second problem is the incorrect use of model parameters. For example, the elastic or viscoelastic modulus of elasticity in a one-dimensional constitutive relation is different from the corresponding elastic or viscoelastic shear modulus in the three-dimensional constitutive equation. Thus, the elastic modulus k in the one-dimensional constitutive relation should be replaced by the elastic shear modulus G and bulk modulus K determined by k and μ in the three-dimensional constitutive relation [40]. (3) When establishing the three-dimensional creep equation, the stress under uniaxial compression is directly replaced by the deviatoric stress S ij under triaxial compression. The relationship between the deviatoric stress S ij and the long-term strength of rock under three-dimensional stress is used as the basis for judging whether the rock enters the third stage [41]. The rationality of this criterion is also questionable. To avoid these problems, a reasonable three-dimensional creep equation is derived below.
In a three-dimensional stress state, assume that the total strain of the nonlinear viscoelastic-plastic creep damage model is ε ij , the strain of the elastomer is ε 1 ij , the strain of the viscous element is ε 2 ij , the strain of the nonlinear Kelvin element is ε 3 ij , and the strain of the viscoelastic-plastic element considering damage is ε 4 ij . As all four elements are in series, the total strain is given as [42]: The stress σ ij at any point in the rock in an elastic state is composed of spherical stress tensor σ m δ ij and deviatoric stress tensor S ij . Similarly, the strain at any point is composed of spherical strain tensor ε ij and deviatoric strain tensor e ij ; hence, we have the following equation [43,44]: ( where δ ij is the Kronecker tensor.
In accordance with the generalized Hooke's law, ( where K is the bulk modulus of the elastomer, MPa; and G 1 is the shear modulus of the elastomer, MPa. Therefore, the three-dimensional constitutive equation of the elastomer is [34]: The three-dimensional creep equation of the viscous element and the nonlinear Kelvin element can be obtained by imitating the creep equation under the one-dimensional stress state: where η 1 ' is the viscous shear coefficient of the viscous element, MPa�h; G 2 is the shear modulus of the nonlinear Kelvin element, Mpa; and η 20 ' is the viscous shear coefficient of the nonlinear Kelvin element, Mpa�h.
The three-dimensional creep relation of the viscoelastic-plastic element considering damage is given by: Where η 30 ' is the viscous shear coefficient of the viscoelastic-plastic element considering damage, Mpa�h; F is the rock yield function; F 0 is the initial value of the yield function (to simplify the calculation, F 0 = 1 [45]); ϕ(�) is a power function, which typically takes the power exponent n = 1 [46]; Q is the plastic potential function, which follows the associated mobility guidelines, namely F = Q; and G 3 is the shear modulus of the viscoelastic-plastic element considering damage, MPa.

Model validation and parameter identification
The rationality and accuracy of the derived model were verified through a rock creep test under freeze-thaw cycles, and the model parameters were determined using mudstone triaxial creep test data. To facilitate the data fitting, the model creep equation was simplified to Eqs. (21) and (22). The fitting results are illustrated in Fig 3, and the nonlinear viscoelastic-plastic creep damage model parameters are stated in Table 1.  The mechanism of rock freeze-thaw damage can be explained as follows. During the freeze-thaw cycle, the freezing of the water inside the rock expands the pores, and the volume expansion generates permanent cracks in the rock as the water melts. Consequently, numerous pores are created inside the rock, damaging the rock and weakening its mechanical properties. Therefore, with the gradual increase in the axial load and the increase in the number of freezethaw cycles, the nonlinear deformation characteristics of rock become more obvious while macroscopic cracks are generated faster until damage occurs. Fig 5 illustrates the quantitative relationship between freeze-thaw times under axial loads and parameters G 1 , G 2 , η 1 ', η 20 ', and λ. The analysis shows that with the increase in the number of freeze-thaw cycles, G 1 , G 2 , and η 20 ' decrease exponentially whereas λ increases exponentially. Hence, the freeze-thaw cycle considerably influences creep deformation, and the influence gradually tends to be stable as the number of freeze-thaw cycles increases. The effect of freeze-thaw times on η 1 ' can be expressed by a cubic polynomial: y = ax 3 + bx 2 + cx + d. Table 2 shows the fitting function of each parameter with changes in the number of freezethaw cycles.

Conclusion
In this study, we connected an elastomer, a viscosity element, a Kelvin element, and a viscoelastic-plastic element in series to construct a nonlinear viscoelastic-plastic creep damage  model. The material parameters and damage factor were also introduced into the model. Oneand three-dimensional creep equations were derived, and triaxial creep data were used to determine the model parameters and to validate the model. The major findings are summarized below: 1. A nonlinear viscoelastic-plastic creep damage model that effectively models the nonlinear creep characteristics of rocks under freeze-thaw cycles was constructed by connecting an elastomer, a viscoelastic element, a Kelvin element, and a viscoelastic-plastic element in series.
2. The established model better reflects the creep behavior of rock in the attenuation creep stage, constant creep stage, and accelerated creep stage. The quantitative relationship between the parameters and the number of freeze-thaw cycles was obtained for the third stage.
3. The nonlinear viscoelastic-plastic creep damage model is simple and reliable, and the parameters are few and easy to control. The model curve and test data have a good fit, with a fitting correlation coefficient close to 1, indicating that the model effectively describes the creep characteristics of rock under different freeze-thaw cycles and verifying the rationality of the model. 4. The quantitative relationship between model parameters and freeze-thaw cycles can be characterized as follows: with the increase in the number of freeze-thaw cycles, parameters G 1 , G 2 , and η 20 ' decrease exponentially, whereas parameter λ increases exponentially. The effect of freeze-thaw times on parameter η 1 ' can be expressed by a cubic polynomial: y = ax 3 + bx 2 + cx + d.